
function [edge_img filtered_img zrc_img] = ecen5283_edge_detect_log(src_img, sigma, hsize, thresh)

%% extract information from src image
src_img = double(src_img);
[rows cols] = size(src_img);

%% Laplacian of Gaussian 
% Generate LoG filter
    function h = gen_log_filter(sigma, hsize)
        K = 2 * hsize + 1;
        h = zeros(K,K);
        for x = 1:K
            for y = 1:K
                u = x - hsize - 1;
                v = y - hsize - 1;
                h(x,y) = (1/(2*sigma^6)) * (u^2 + v^2 - 2*sigma^2) * exp(-(u^2 + v^2)/(2*sigma^2));
            end
        end
%         h = h - mean(h(:));
%         h = h * sigma^2;
    end
% sigma = 2;
% hsize = ceil(sigma * 3);
hlog = gen_log_filter(sigma, hsize);
% hlog = fspecial('log', 2*hsize+1, sigma);
% LoG = [0 0 1 0 0;0 1 2 1 0;1 2 -16 2 1; 0 1 2 1 0; 0 0 1 0 0];

% Filter image by LoG filter
filtered_img = imfilter(src_img, hlog, 'conv');
zc_img = zeros(rows, cols);
edge_img = zc_img;
%% Find zero-crossing
for i = 2:rows-1
    for j = 2:cols-1        
         % process boundary - nearest padding
         up_i = max(1,i-1);
         down_i = min(rows,i+1);
         left_j = max(1,j-1);
         right_j = min(cols,j+1);

         % compute 8 neighbor pixels
         n = filtered_img(up_i,j);             
         s = filtered_img(down_i,j);
         w = filtered_img(i,left_j);
         e = filtered_img(i,right_j);
         nw = filtered_img(up_i,left_j);
         ne = filtered_img(up_i,right_j);
         sw = filtered_img(down_i,left_j);
         se = filtered_img(down_i,right_j);

         %zero-crossing check                  
         max_zc = 0;
         if (n*s < 0)
             max_zc = abs(n) + abs(s); 
         elseif (e*w < 0)
             max_zc = max(max_zc,abs(e) + abs(w)); 
         elseif (ne*sw < 0)
             max_zc = max(max_zc,abs(ne) + abs(sw)); 
         elseif (nw*se < 0)
             max_zc = max(max_zc,abs(nw) + abs(se));          
         end
         zc_img(i,j) = max_zc;
    end
end
%% Thresholding
edge_img(zc_img >= thresh) = 1;
edge_img(zc_img <  thresh) = 0;

%% Post processing
edge_img = im2bw(edge_img);
zrc_img = edge_img;
edge_img = remove_isolated_point(edge_img);
edge_img = bwmorph(edge_img,'thin');
end